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ABSTRACT 



^ Aims. A complete set of diagnostic tools aimed at producing synthetic synchrotron emissivity, polarization, and spectral index maps from 
relativistic MHD simulations is presented. As a first application we consider here the case of the emission from Pulsar Wind Nebulae (PWNe). 
OO Methods. The proposed method is based on the addition, on top of the basic set of MHD equations, of an extra equation describing the evolution 
1 of the maximum energy of the emitting particles. This equation takes into account adiabatic and synchrotron losses along streamlines for the 
distribution of emitting particles and its formulation is such that it is easily implemented in any numerical scheme for relativistic MHD. 
Results. Application to the axisymmetric simulations of PWNe, analogous to those described by Del Zanna et al. (2004, A&A, 421, 1063), 
allows direct comparison between the numerical results and observations of the inner structure of the Crab Nebula, and similar objects, in the 
optical and X-ray bands. We are able to match most of the observed features typical of PWNe, like the equatorial torus and the polar jets, with 
^ ■ velocities in the correct range, as well as finer emission details, like arcs, rings and the bright knot, that turn out to arise mainly from Doppler 
CIh" boosting effects. Spectral properties appear to be well reproduced too: detailed spectral index maps are produced for the first time and show 
q . softening towards the PWN outer borders, whereas spectral breaks appear in integrated spectra. The emission details are found to strongly 
^ 1 depend on both the average wind magnetization, here cr eff & 0.02, and on the magnetic field shape. 

0^ | Conclusions. Our method, in spite of its simplicity, provides a realistic modeling of synchrotron emission properties, and two-dimensional 
Ci ■ axisymmetric relativistic MHD simulations appear to be well suited to explain the main observational features of PWNe. 
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1- Introduction loud quasars, to the bulge and disk of our Galaxy, and galac- 
tic sources like Supernova Remnants (SNRs). Precisely to a 

Non-thermal radiation, and synchrotron emission from rela- subclass of the latter belongs one of the most studied astro _ 

tivistic particles in magnetized plasmas in particular, is a ubiq- physical sources: the Crab Nebula . The Crab Nebula is the 

uitous signature in high energy astrophysics. Strong shocks in prototype of the S0 - C alled Pulsar Wind Nebulae (PWNe), or 

relativistic flows are believed to be the natural site for the ac- plerionS9 and it is the object where synchrotron radiation was 

celeration of particles up to relativistic velocities. High energy first identified (Shklovsky E953J. The emission from PWNe 

electrons and positrons, then, efficiently emit synchrotron ra- displays all the characteristic properties of synchrotron radi- 

diation while spiraling around magnetic field lines. For opti- ation from re i at ivistic particles: a continuous, very broad-band 

cally thin sources, the emission spectrum, typically a power spectrum> extending from radio to X-rays and beyond (though 

law, is related to the particles' distribution function at injection Inverse Compton typica n y dominates at gamma-ray frequen- 

through well known transport and radiation physics. However, des) and a high degree of lmear po i arization; spec tral indices, 

the connection may be not easy to establish quantitatively in the defined by Fy K y -a where Fy is the net flux? range between 

presence of multi-dimensional flow patterns, making it difficult < a < 1 .2 and in a given source usually steepen with increas- 



ing frequency, 



to extract from observations direct information on the most ba- 
sic physics of the sources. 

Environments where synchrotron emission is encountered While the theoretical modeling of the emitted radiation 

span from extragalactic sources, like cluster halos or Active soon reached a high level of sophistication (e.g. Kennel & 

Galactic Nuclei (AGNs), for example jets and lobes in radio- Coroniti 1984b for PWNe), what has been lacking until very 

recently is an appropriate modeling of the fluid and magnetic 

Send offprint requests to: Luca Del Zanna structure of the relativistic plasma responsible for the emis- 

e-mail: ldz@arcetri.astro.it sion. In the last few years, however, there has been an in- 
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creasingly rapid growth in accuracy and robustness of shock- 
capturing numerical schemes for relativistic hydrodynamics 
and MHD (see Marti & Muller l2003l and Font 120031 for peri- 
odically updated reviews), so that nowadays relativistic simu- 
lations are becoming a standard investigation tool for plasma 
astrophysics. Among the others, we mention here the works on 
special relativistic MHD by Komissarov ( 1999 ) and Del Zanna 
et al. (2003), with the latter putting forward a simpler and yet 
powerful alternative to characteristics-based algorithms. Based 
on these techniques, conservative shock-capturing schemes for 
MHD in General Relativity (GRMHD) have also been con- 
structed (e.g. Gammie et al. 2003 ; Komissarov 2005 ; Anton et 
al. 120061) . especially for studies of the accretion flows around 
Kerr black holes or extraction of gravitational energy and sub- 
sequent jet acceleration. Recently, coupling of conservative 
schemes for GRMHD with solvers for the Einstein equations 
has been achieved too (Duez et al. 2005 ), for studies of binary 
neutron star mergers, core collapse in supernovae, black hole 
formation and gamma-ray burst generation. 

The goal of the present paper is to provide the numeri- 
cal community with a simple computational tool that can be 
used to compute the brightness, polarization and spectrum of 
the synchrotron light associated with the plasma dynamics re- 
sulting from a given MHD simulation. The method proposed 
is quite general and can be applied to any scheme for relativis- 
tic MHD (not necessarily in conservative form, either in full 
3-D or under particular symmetries) and it just requires the nu- 
merical solution of an additional evolution equation, as we will 
discuss later. As a first application of our techniques we study 
here in detail the emission properties of PWNe, as they can 
be computed based on the recent axisymmetric model by Del 
Zanna et al. (120041 . 

The interest of the scientific community for PWNe has re- 
cently received considerable impulse from the X-ray observa- 
tions made by Chandra and XMM-Newton, showing that the 
Crab Nebula and other PWNe all present complex but similar 
inner structures: a higher emission torus in what is thought to be 
the equatorial plane of pulsar rotation, typically with brighter 
features like arcs or rings, a central knot, and one or two oppo- 
site jets along the polar axis, apparently originating from very 
close to the pulsar position (Hester et al. 119951 for Hubble and 
ROSAT earlier observations; Weisskopf et al. 2000; Hester et 
al. 120021 Helfand et al. I200T1 Pavlov e t al. I200T1 Go tthelf & 
Wan gl2000l Gaensler et al.l200Tll2002l Lu et al.l2002l Camilo 
et al. 120041 Slane et al. 120041 Romani et al. 120051 . This sce- 
nario is generally referred to as the jet-torus structure. 

While the presence of an equatorial torus of enhanced 
emission can be explained by applying the standard 1-D rel- 
ativistic MHD solutions (Kennel & Coroniti 1984a Emmering 
& Chevalier 1987) and assuming that the pulsar wind has 
a higher equatorial energy flux (Bogovalov & Khangoulian 
2002), the origin of the polar jets has puzzled plasma physi- 
cists until recently, given the known difficulties in obtaining 
self-collimation via hoop stresses by the toroidal field in an 
ultra-relativistic MHD outflow (Lyubarsky & Eichler 2001 and 
references therein). However, an axisymmetric energy flux en- 
hanced at the equator produces an oblate termination shock 
with a cusp-like shape at the poles. Moreover, if hoop stresses 



are at work in the mildly relativistic post-shock flow, jet col- 
limation may occur on the axis at these cusps, thus giving 
the impression of an origin inside the wind itself (Lyubarsky 
2002 ). This scenario has been recently confirmed by relativistic 
MHD numerical simulations (Komissarov & Lvubarskv 120031 
120041 Del Zanna et al. 12004ft . which all show that the ter- 
mination shock really assumes the predicted shape, polar jets 
with velocities in the observed range 0.5 - 0.8 c are formed 
(provided the wind magnetization is high enough (approxi- 
mately 1%), and the complex flow patterns which arise near 
the termination shock are responsible for most of the observed 
brighter features, basically just due to Doppler boosting effects 
(Komissarov & Lyubarsky 2004). 

Compared to the quality of the simulation results for the 
flow structure, synchrotron emission modeling has been so far 
rather crude. In the only attempt so far at computing syn- 
chrotron surface brightness maps on top of the flow struc- 
ture obtained from 2-D MHD simulations (Komissarov & 
Lyubarsky 2004 ) radiative losses were approximated by an ex- 
ponentially decaying factor, with a spatially constant cooling 
time depending on the frequency of observation and an aver- 
age magnetic field. A power law distribution function for the 
emitting particles, normalized to the local thermal pressure to 
take into account variation of specific volume, was assumed in 
the emissivity. A more sophisticated method was proposed by 
Shibata et al. (2003 ): here several particle energies are evolved 
independently in order to reconstruct the initial distribution 
function at any location and time. However, this method is suit- 
able only for nearly steady- state solutions where the individual 
paths (the fluid streamlines) can be unambiguously followed. 
The method proposed here is a step forward concerning the ap- 
plication to time-dependent numerical schemes: emitting par- 
ticles are continuously injected at the termination shock and 
then the full particles' energy equation is integrated in time. 
However, differently from the the above cited paper, just the 
maximum energy attainable is evolved, and it is taken as the 
upper cut-off in the local power law distribution function that 
enters the emission coefficient. 

Other than just synchrotron surface brightness, we also 
show here how to build polarization maps, namely maps of 
Stokes parameters, degree of linear polarization, and position 
angle. Such data are available for the Crab Nebula in radio 
(Wilson [T9721 Velusamy [T985l and optical (Schmidt et al. 
[T9791 Hickson & van der Berg IT9901 Michel et al. [T99ll 
though high resolution data are still missing. X-ray polariza- 
tion maps are not available yet. Synthetic maps have been re- 
ported in a preliminary work (Bucciantini et al. 2005b), where 
optical maps have been produced for one test case, here we 
discuss the dependence on the model parameters. Finally, our 
method allows us to produce detailed spectral maps, which 
can be compared to optical observations (e.g. Veron-Cetty & 
Woltier fT993l and to Chandra and XMM-Newton X-ray data 
(Weisskopf et al. l2000l Willingale et al. l200Tl Mori et al.[2004 ). 
Synthetic spectral index maps can provide us with an additional 
diagnostic tool for the inner magnetic structure of PWNe, since 
the spectral softening due to synchrotron losses depends on the 
local magnetic field. To our knowledge, this is the first time 
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synthetic synchrotron spectral index maps based on data from 
numerical simulations have been produced. 

The paper is structured as follows: the proposed method is 
described in Sect. |3 the main PWN model assumptions are 
summarized in Sect. |3 numerical results and synthetic syn- 
chrotron maps are shown in Sect. |4l where also a comparison 
with observations is attempted. Finally, Sect. is dedicated to 
the conclusions. 



2. Synthetic synchrotron emission recipes 

In the present section a method for producing synthetic syn- 
chrotron surface brightness and polarization maps, based on 
the results of numerical simulations, will be derived and dis- 
cussed. The proposed method takes into account both adiabatic 
and synchrotron losses in the evolution of the emitting particles 
advected by the nebular flow, and it can be easily implemented 
in any numerical scheme for relativistic MHD equations. The 
only assumption made here is that the code provides the vari- 
ables p (proper mass density), v (flow bulk velocity), p (proper 
thermal pressure), and B (magnetic field) as functions of space, 
with or without special symmetries reducing the dimension- 
ality of the system, and time. Lorentz transformations lead- 
ing to relativistic corrections in the emission properties, like 
Doppler shift and boosting or polarization angle swing, will 
also be fully considered. Therefore, in the following subsec- 
tions, different notations will be employed where ambiguities 
may occur: primed quantities will refer to the reference frame 
comoving with the fluid, whereas standard notation will indi- 
cate quantities in the observer's frame. 

2.1. Distribution function evolution 

The observed spectra in astrophysical sources of non-thermal 
radiation can be often approximated by power laws, with vary- 
ing spectral index a in different frequency regimes. It is there- 
fore natural to model the emitting particles' (typically electrons 
and positrons) distribution function at injection locations as a 
combination of power laws as well. Let then 6o be the particle's 
energy (in units of the rest mass energy mc 2 ) and consider a dis- 
tribution function which is isotropic in momentum and power 
law in energy, within a given energy range: 



and from the above definitions we also find 



r ( \ ^ -(2a+l) 



< € < e£ 



(1) 



where from now on quantities labeled with the index will 
refer to the injection region. The constant of proportionality 



Po 



A = K n n = K p — - 



(2) 



can be determined by the definitions of the proper number den- 
sity no and thermal pressure po in terms of fo(eo), provided 
e™ ax » 6™ in (see Kennel & Coroniti [T984bl) . If a > 1/2, as it 
is the case in the following, integration over momentum space 
and energy gives 



K n = 2a{e^) 2 \ 

K p = 3(2a-l)(e™») 2 <*-\ 



(3) 



e o 



2a- 1 



Po 



2a / nomc 2 



(4) 



that can be determined by the Rankine-Hugoniot relations if the 
injection location is a shock (6™ m is of order of the upstream 
wind Lorentz factor, ~ 10 6 for PWNe). 

Time evolution of the single energies along post-shock 
streamlines is governed by adiabatic and radiative losses. If 
d/df = y(d/dt + v • V) is the total time derivative in the co- 
moving frame, the equation for e is 



Ai n6= Aw/3 

df df 



1 

+ - 

6 



d_6 

df 



sync 



where synchrotron losses are 

de\ 4e 4 
df 



sync 



9m? c 



(5) 



(6) 



in which we have averaged over all possible pitch angles with 
the local magnetic field that a particle may experience during 
its lifetime. Based on the evolution of the single energies, we 
now seek an equation for the evolution of the distribution func- 
tion /(e) along streamlines. This problem was solved, for a sta- 
tionary radial flow, by Kennel & Coroniti ( 1984b), extended to 
the multi-dimensional case by Begelman & Li ( 1992), and first 
applied to the results of a time-dependent code by Bucciantini 
et al. (12005 al) . From the integration of Eq. (|5|) along streamlines 
it is possible to write 



J/3 



n l/3 
n o 



,1/3 



(7) 



where 6co takes into account the integrated synchrotron losses 
and corresponds to the remaining energy of a particle with 
energy 6o —> oo at the injection location, that is to say the 
local maximum energy attainable. By using the above rela- 
tion it is possible to express 6 as a function of eo and 6oo, 
which can be easily calculated for a stationary flow every- 
where in space. From the conservation of the particles' number 
f(e)de/n - /o(£b)deoMo along streamlines we are finally able 
to derive the distribution function as 

n 2a-l 



P P c -(2a+l) 



4n mc z 



1/4 



,6 < 6oc 



(8) 



where we have used adiabaticity of the relativistically hot 
plasma: p ~ n 4/3 . 

The method we propose in the time-dependent case, pro- 
vided changes in the global structure occur on timescales not 
shorter than the typical transport times, is to evolve also the 
variable 6oo by adding an extra equation to the relativistic MHD 
code. For a shock-capturing scheme the required conservative 
form of Eq. © is found by combining it with the continuity 
equation for p = nm: 



(9) 



|( r p£) + V. (r p£ D ) = ^fe) , 

dt eoo \ df / sync 

with £ = 6oo/p 1/3 and e M being a large enough number at in- 
jection, so as to mimic the theoretical behavior over the largest 
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possible domain (differences due to €q < oo should be confined 
to the close vicinities of the injection location, see Kennel & 
Coroniti 1984bJ. 

Once we have £00 as a function of time and space, we can 
model the effects of adiabatic and synchrotron losses on the 
distribution function by applying Eq. ®. For ease of imple- 
mentation in a numerical code, however, a further approxima- 
tion is finally needed. Since no and po are difficult to trace 
back starting from the position along the local streamline, in 
the multidimensional case one should also assume that their in- 
fluence is negligible. This occurs when two conditions are met: 
the dependency on the term (p/po) l/4 is weak and so is that on 
quantities involving the temperature ~ poM) at injection. Both 
conditions are easily satisfied for 2a - 1 ^ 0, so that K n and K p 
are approximately constant. 

2.2. The emission coefficient 

Consider an ultra-relativistic electron (or positron) with nor- 
malized energy 6, which spirals around the local magnetic field 
B f . Its synchrotron spectral power (per unit frequency) is given 
by 



2e 4 
3m 2 c 3 



5lVS(v',0, 



(10) 



where B' ± is the field component normal to the particle's veloc- 
ity, the spectral density is 



s{vy c ) 



8nv' Ivi 



(11) 



and the function F(x) (see Rybicki & Lightman 1979 for its 
properties) is such that emission peaks around the critical fre- 
quency 



3e 
4nmc 



B\e 2 . 



(12) 



The local emission coefficient f v (radiated power per unit fre- 
quency, volume and solid angle), for given observation fre- 
quency V and direction n', is obtained by integrating the prod- 
uct of Eqs. <[T0]> and ® over the whole range of local particle 
energies: 



f v (v',n') =j P(/,e)f(e)de. 



(13) 



Note that, since radiation from ultra-relativistic particles is 
strongly beamed in their instantaneous direction of motion 
(within a cone of angle 1/6 <?c 1), only those particles with 
pitch angle coinciding with the angle between B' and n' (thus 
= \B' x n'\) contribute to the emission along the line of sight 
in Eq. O- 

The above integral can be calculated by replacing the en- 
ergy with the variable x - V jV c oc e~ 2 . In order to take advan- 
tage of the known properties of the function F(x) we further 
simplify Eq. © to a standard power law form 



/(6) = _£^ 6 -(2-i), 
47T mc L 



e < 60c 



(14) 



in which we have assumed that K p is constant, the dependency 
on the term (p/po) l/4 is weak, and synchrotron burn-off can be 
modeled with a sharp cut-off at 600 . This last additional assump- 
tion is again a good approximation provided 2a - 1 ^ 0. The 
result is 



f v (v',ri) = Cp\B'xn' 



,f\a+lyf-a 



(15) 



where all spatially independent terms have been gathered in the 
constant 



c _ V3 a+5/3 / of+S/3 \ r / a + 1 /3 
167T a + 1 



3e 



m 2 c 4 \ 27imc \ 



K p , (16) 



in which the property T(£ + 1) = ^r(^) of the Euler gamma 
function T(£) has been used. 

In order to obtain the emissivity in the observer's fixed 
frame of reference, relativistic corrections must be taken into 
account. It is well known that both the frequency and the emis- 
sion coefficient itself transform via the Doppler boosting factor 



D = 



1 



y(l-j8-n) 



(17) 



as v = DV and j v = D 2 f v , where f$ = v/c. The resulting emis- 
sion coefficient is finally written as 



Cp \B' x n'\ a+1 D a+2 v- a , Voo > y, 



7v(v, n) 



(18) 



(19) 



(0, Voo < y. 

The cut-off frequency for synchrotron burn-off is 

Voo = DV^) = D-p- \B' x n'\ei, 

where Doppler shift has been considered and where 600 is 
provided directly by the numerical scheme via integration of 
Eq.Q. 

The terms B' in Eq. (|6|), needed to evolve £00, and \B' x n'\ 
in Eqs. (1181191) . are defined through quantities in the comov- 
ing frame. Application of the composition rule for relativistic 
velocities to the observer direction versor n yields 



D 



n + 



7+1 



(20) 



while Lorentz transformations of Maxwell's equations in the 
(ideal) MHD case E + fi x B = provide 



B f 



B + 



r 



7 + 1 



(P-B)fi 



(21) 



After some straightforward algebra, the required quantities can 
be written as 



B' = - Jb 2 +j 2 (P-B) 2 , 



(22) 



and 



\B f x n'\ = - Jb 2 - D\B • n) 2 + 2yD(B • n)(J3 • B\ (23) 
r v 

in which everything on the right-hand- side is now measured in 
the observer's frame. 
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2.3. Surface brightness and polarization maps 

To compare the synthetic synchrotron emission obtained by nu- 
merical simulations with real images, the first step certainly 
requires the production of surface brightness maps. Let us con- 
sider a Cartesian observer's reference frame in which X lies 
along the line of sight n and Y and Z are in the plane of the sky 
(say with Z indicating the North). By switching from the lo- 
cal coordinate system, which may present special symmetries, 
to the observer's system (X, 7, Z) (taking into account the an- 
gles of inclination of the emitting object), surface brightness 
(or specific intensity) maps are obtained as 



/v(v,F,Z) 



-r 



j v (v,X,Y,Z)dX. 



(24) 



A powerful additional diagnostic tool may be provided 
by polarization measurements, where available. Synchrotron 
emission from relativistic particles is known to be linearly po- 
larized, with a high degree of polarization, and Stokes param- 
eters U v and Q v (V v = 0) can be calculated from numerical 
simulations by measuring the local polarization position angle 
X, that is the angle of the emitted electric field vector e in the 
plane of the sky, measured clockwise from the Z axis. Synthetic 
maps of the Stokes parameters are thus built as 

Q v (y,Y,Z) = — j y (v,X,Y,Z)cos2 X dX, (25) 

a + 5/3 J_ TC 

ex ~\~ i r* 00 

U v (y,Y 9 Z)= — Mv,X 9 Y 9 Z)sm2xdX, (26) 

a + 5/3 J.oo 

where the additional factor comes from the intrinsic proper- 
ties of synchrotron emission (see Rybicki & Lightman 1979). 
Related to the Stokes parameters above is the degree of (linear) 
polarization, or polarization fraction, defined as 



n v 



ly 



(27) 



In order to calculate^, relativistic effects like position angle 
swing must be taken into account (Blandford & Konigl fT979l 
Bjornsson 1982 Lyutikov et al. 2003 ). In the comoving frame, 
the emitted electric field of the linearly polarized radiation is 
normal to both the local magnetic field B f and the line of sight 
n\ thus e' oc n'xB', and the radiated magnetic field is therefore 
b' - n' x e' . Lorentz transformations give the electric field in 
the observer frame as 



r 



y+l 



(e'J3)J3-J3xb' 



(28) 



and by using Eqs. (f20b and (THI) we can express the required 
electromagnetic fields in terms of quantities defined in the ob- 
server frame. After some calculations, the simple formula de- 
rived by Lyutikov et al. ( 2003 ) is retrieved: 



e oc nxq, q = B + nxifixB). 



(29) 



Notice that e lies in the plane of the sky, as expected, and that 
q correctly reduces to B for non-relativistic fluid bulk motions. 



We are now ready to derive the position angles in terms of the 
components of q in the plane of the sky, which are 

qy = (1 -Px)By +PyB x , qz = (1 ~Px)B z + PzB x . (30) 

The functions of the position angle appearing in Eqs. (1251261) 
may be written in terms of these as: 



q\-q 2 z . 2q Y qz 
cos 2x= — J' sm2 * = — 2 2* 



<t Y + <tz 



<fi + <tz 



(31) 



2.4. Spectral index maps and integrated spectra 

In order to complete our set of diagnostic tools based on the 
synthetic synchrotron emission obtained from numerical simu- 
lations, we now turn to spectral properties. An interesting pos- 
sibility is provided by maps of spectral index. This quantity is 
defined as 



or v (vi, V2, Y,Z) = - 



log[I v (v 2 , Y,Z)/I v (v u Y,Z)] 
log(v 2 /vi) 



(32) 



for any given interval in frequency [vi, V2L In general, due to 
synchrotron radiative losses, this quantity is expected to in- 
crease, thus the spectrum is expected to steepen, when mov- 
ing away from the source of ultra-relativistic particles (spectral 
softening). 

Finally, integrated spectra (or net fluxes) are obtained by 
averaging the specific intensity over the solid angle covered by 
the source in the plane of the sky, which may be calculated as 



F v (v) 



Y,Z) d7dZ, 



(33) 



where d is the distance of the astronomical object under con- 
sideration. 

Notice that both the above quantities related to spectral fea- 
tures have a meaning only if synchrotron burn-off is modeled 
in some way by the numerical simulation. 

3. Pulsar wind model and initial settings 

As anticipated in the Introduction, let us now apply our set 
of diagnostic tools to the modeling of the inner structure of 
PWNe. Axisymmetric 2-D simulations of the interaction of an 
ultra-relativistic magnetized pulsar wind with supernova ejecta 
expanding in the static interstellar medium were described in 
Del Zanna et al. (2004), to which the reader is referred for all 
the modeling details and for the results. Here we will calcu- 
late the synchrotron emission from the results of new simu- 
lations, obtained by using very similar settings. The numeri- 
cal code employed is described in Del Zanna & Bucciantini 
(l200a , Del Zanna et al. (120031) . and Londrillo & Del Zanna 
(2004). It is a shock-capturing code solving the ideal relativis- 
tic MHD equations in conservative form, to which Eq. © has 
been added. Note that the effects of synchrotron radiative losses 
on the fluid total energy evolution are not taken into account. 
Radiative losses would become important if the particles radi- 
ated away a considerable fraction of the pulsar energy input to 
the nebula. However, the radiated power is estimated to be of 
order 10% of the instantaneous pulsar input in the case of the 
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Crab Nebula, and generally much smaller than this in the case 
of other PWNe, so the condition of adiabaticity should gener- 
ally be considered as a good approximation. 

Let us briefly summarize the initial conditions and the pul- 
sar wind model in particular. The numerical box is made up of 
400 cells in the radial direction, between r = 0.1 ly and r = 10 
ly (a logarithmic grid with 200 cells per decade is used to en- 
hance resolution in the inner region), and 100 cells in the latitu- 
dinal direction, with polar angle 6 defined in the first quadrant 
only (symmetry conditions are applied at the equator 6 = n/2). 
Within an arbitrary radius of 0.2 ly we impose initially the wind 
conditions, while the expanding ejecta, with linearly increas- 
ing radial velocity, occupy the region up to 1 ly. The (constant) 
wind luminosity is taken to be Lq = 5 x 10 38 erg s" 1 , the en- 
ergy of the supernova explosion turned into kinetic energy of 
the ejecta E Q j = 10 51 erg, and the total mass of the ejected ma- 
terial M e j = 3M©. For similar settings and a description on the 
(1-D) evolution of the PWN / SNR system see van der Swaluw 
et al. (2001 ) and Bucciantini et al. (2003). The chosen values 
were seen to provide an overall evolution compatible with that 
of the Crab Nebula. Furthermore, the magnetic field B is as- 
sumed to retain just the toroidal component, as in Kennel & 
Coroniti ( 1984a), while the velocity v is assumed to be always 
poloidal, so that some of the relativistic corrections in Sect. 12.21 
are simplified since fi • B = 0. 

The latitudinal dependence of the wind energy flux is pro- 
vided by the following choice for the wind bulk Lorentz factor 



0.05 p 



wind magnetization 



y(6) = yo [<*o + (1 - <*o) sin 2 o] , 



(34) 



where 70 = 100 is the equatorial value and = 0.1 the 
anisotropy parameter. Velocity is assumed to be in the radial 
direction alone at t = 0. The residual magnetic field at large dis- 
tances from the pulsar's light cylinder is thought to be mainly 
toroidal with a dependence B ~ sinO/r (split monopole mod- 
els: Michel [T^731 ) inside the wind. By prescribing such a field 
and an isotropic mass flux, the energy flux for a cold ultra- 
relativistic wind may be then written as 



F(r, 6) = F J [ a + (1 - a + &) sin 2 0] , 



(35) 



where cr is the wind magnetization parameter, namely the ratio 
of magnetic to kinetic energy fluxes computed at the equator 
(Kennel & Coroniti 1984a), and Fo is the equatorial kinetic 
energy flux at the reference radius ro, to be derived in terms of 
L , c*o, and cr by integrating Eq. (l35l over a spherical surface 
of radius r. 

The value of cr is typically expected to be a small number 
in order to match the predicted and observed nebular emission 
(Kennel & Coroniti [T^84bl cr = 3 x 10" 3 for the Crab Nebula), 
thus one is left with the problem of how to convert a Poynting 
flux dominated outflow near the pulsar to a kinetically domi- 
nated wind before the termination shock (the so-called sigma 
paradox) (see e.g. Arons 2004). One of the solutions proposed 
invokes the presence a striped wind region around the equa- 
tor, due to the inclination of the pulsar's magnetosphere with 
respect to the rotation axis, where dissipation in the resulting 
current sheets may form a neutral (unmagnetized) equatorial 
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Fig. 1. The wind magnetization B 2 /4npc 2 y 2 as a function of 
the polar angle 6, at any radius r, for the two runs defined in 
Eq. 03- In both cases the averaged value is approximately the 
same: <x e ff ~ 0.02. 

region (Coroniti 1990). We model here such a situation by pre- 
scribing the wind toroidal field as 



B(r, 



sin 6 tanh 



.(H 



(36) 



where we set Bo = ^4ncrFo/c and where the last term mim- 
ics the conversion of magnetic to kinetic energy in the striped 
wind. The parameter b controls the angular size of the zone of 
lower magnetization: a large b indicates a narrow striped wind 
region, and the ideal limit is recovered by letting b — > oo. Here 
we assume that Eq. (l35l still holds even for a finite b. Finally, 
the density p, decreasing as r~ 2 , can be derived from the above 
equations and will depend on cr and b too, since it is supposed 
to account also for the conversion of magnetic to kinetic energy 
in the striped wind region (note that the mass flux is actually 
isotropic only in the ideal limit b — > oo). 

As our main goal here is the comparison with observations 
of the Crab Nebula, the temporal evolution of the PWN / SNR 
interaction will not be followed in any detail. On the contrary, 
we will focus our attention on the synchrotron emission proper- 
ties of the simulated PWN at a time corresponding to the age of 
the Crab, t « 1000 yr, which is still in the free expansion phase 
(Bucciantini et al. 120031) . As shown by Del Zanna et al. ( 2004 ), 
equatorial flows and polar jets with velocities in the range of 
the observed values are found for an effective magnetization of 
the wind plasma <x e ff ^ 0.01, where <x e ff is obtained by aver- 
aging over polar angle (the magnetization in the wind region 
does not depend on r). Thus, here we will consider two differ- 
ent runs, A and B, with approximately the same <x e ff ~ 0.02 
and very different magnetic fields: 



run A : cr = 0.025, b = 10, 
run B : cr = 0.1, b = 1. 



(37) 



The wind magnetization B 2 /4npc 2 y 2 for the two runs as a func- 
tion of the polar angle 6 is shown in Fig. [T] Notice the differ- 
ence in the striped wind region width around the equator: run A 
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refers to a nearly ideal case, as those discussed by Del Zanna et 
al. (|2004), while in run B the striped wind region is very wide. 

4. Simulation results 

4.1. Flow structure 

Before studying the synchrotron radiation maps, let us briefly 
comment on the flow structure and related properties for the 
two different runs, focusing on the quantities that are relevant 
for computing the emission. In Fig. |2| we show the velocity 
field (velocity magnitude and streamlines), the magnetization 
B 2 /4nwy 2 (w = pc 2 + 4p is the enthalpy, dominated by the 
kinetic component in the pulsar wind and by the thermal one 
inside the nebula), and the maximum particle energy 6oo. The 
latter is evolved according to Eq. (|9|> and its value at injection 
is 10 9 . Initialization of 6co occurs, like for all other quantities, in 
the wind region, so that the value of 10 9 is further enhanced at 
the crossing of the termination shock, due to the conservation 
of the quantity 8 = 6co/p 1/3 rather than €oo itself. 

The first thing to notice in Fig. |2| is that the wind zone, 
characterized by the central white region with v « c and oblate 
shape, is smaller in run A, because of the stronger pinching 
forces around the equator. This difference in the size of the ter- 
mination shock in the two cases will have visible consequences 
in the emission maps. The flow structure is similar: in both 
cases substantial velocities are present in the equatorial plane 
as well as along the polar axis (the supersonic jets). It should 
be emphasized that the values we find, v « 0.5-0.8 c, are in 
agreement with those observed in the Crab Nebula (Hester et 
al. l2Q02b and in Vela (Pavlov et al. 2003). The complex flow 
pattern around the termination shock is also similar, but hoop 
stresses, which divert the flow toward the polar axis through 
small scale vortices, occur closer to the termination shock it- 
self in case A, due to a much larger region of increased magne- 
tization around the equator. This is clearly visible in the central 
panels: equipartition is reached or even exceeded in a large sec- 
tor around the equator for run A, while this is the case basically 
just along the termination shock front for run B, due to local 
compressions of the magnetic field by the flow vortices. Notice 
that the situation described here is very different from that en- 
countered in standard stationary or self- similar radial models, 
where the nebular plasma reaches equipartition only at large 
radii. This is entirely due to the relaxation of the 1-D approxi- 
mation, and to the complex flow pattern around the termination 
shock determined by the anisotropic wind energy flux. 

The two panels on the right hand side show 6co, that is the 
highest possible energy of emitting particles at a given loca- 
tion. As we can see, in both cases the pattern resembles that of 
velocity, as could be expected given that the particle energy is 
advected by the flow. On the other hand, these rightmost pan- 
els are the "negative" of the corresponding central ones repre- 
senting the magnetization structure. This is due to synchrotron 
losses. And again these losses determine the difference between 
the maps of 6oo for the two runs: while in case A the external 
parts of the map are dark (a lower cut at 6 = 10 6 has been 
imposed), in run B the same regions have a higher maximum 
energy. This different behavior arises from the fact that in run 



A particles loose most of their energy in the high magnetiza- 
tion equatorial channel, whereas in case B, with a lower mag- 
netization there, higher energy particles survive and are then 
transported across the nebula by the large scale vortices. The 
structure of the magnetic and velocity fields and the resulting 
distribution of 6co will have important consequences, especially 
for the spectral properties of synchrotron emission, as will be 
widely discussed in Sect. 14.41 

4.2. Surface brightness maps 

In order to obtain synchrotron emission maps, the spectral in- 
dex a of the power law in Eq. {0 must be assigned. It is well 
known that in the case of the Crab Nebula more than a single 
power law at injection is needed in order to explain the ob- 
served spectrum from the radio to the X-ray band (e.g. Amato 
et al., 2000). Here we will assume that particles responsible 
for optical and X-ray emission come from the same distribu- 
tion function, injected at the termination shock front, while the 
radio emitting particles are not considered. Therefore we take 
a = 0.6, in agreement with the values reported in optical spec- 
tral index maps (Veron-Cetty & Woltjer, 1993 ) close to the ter- 
mination shock, where synchrotron burn-off has not occurred 
yet. Note that this value satisfies the condition 2a - 1 ^ 0, 
needed to approximate the distribution function in Eq. © with 
the simple power law in Eq. (THl . We have verified, by measur- 
ing the ratios p/po and po/no, that the combined effect of the 
neglected terms leads to changes in the emissivity of ~ 10%, at 
most. 

In Fig. surface brightness maps are shown in optical and 
X-rays for runs A and B (in logarithmic scale and normalized 
to the respective maximum value). Optical images are calcu- 
lated for A = 5364 A (left panels), one of the wavelengths 
selected by Veron-Cetty & Woltjer ( 1993), at which emission 
from the outer filaments is negligible. For the X-rays, instead, 
we use hv = 1 keV, a value in the range of both the satellites 
Chandra and XMM-Newton. The images are obtained by ap- 
plying Eq. (l24l and assuming an inclination of the symmetry 
axis of 30° with respect to the plane of the sky and of 48° with 
respect to the North, values appropriate for the Crab Nebula 
(e.g. Weisskopf 2000). A square box of 6 ly x 6 ly, centered on 
the pulsar position, is used for our synthetic images. 

The first thing to notice is that the emission is more dis- 
tributed at optical frequencies compared with the X-rays, where 
the external regions appear darker. This is the expected ef- 
fect of synchrotron burn-off, that causes the sources to appear 
smaller with increasing observation frequency. The central re- 
gions are instead rather similar in the two bands, with evi- 
dence for brighter features like rings, arcs and a central knot. 
The latter, in particular, tends to dominate the overall emission 
(especially in run B) if a linear scale is employed for the im- 
age. As it was first shown by Komissarov & Lyubarsky ( 2004), 
these bright features are the result of Doppler boosted emission 
from fluid elements moving at relativistic speeds toward the ob- 
server. From a comparison with Fig.|2|it is clear that the regions 
with high velocity and magnetization, and thus higher emissiv- 
ity, are those around the termination shock: the multiple rings 
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Fig. 2. Velocity magnitude (in units of c) and streamlines (left panels), magnetization (center panels), and maximum particle 
energy 6co (in units of m e c 2 , right panels) in the inner region of the simulated PWN, for run A (upper row) and run B (lower row). 
Distances from the central pulsar are reported on the axes, expresses in light year (ly) units. 



observed are due to the external vortices produced by the hoop 
stresses which divert the flow, the brighter arc is due to flow 
along the termination shock, while the knot is due to the fluid 
escaping from the polar cusp-like region toward the observer 
(see the above cited paper for a graphical explanation). 

Morphological differences arise from the comparison of run 
A and B, in the X-ray maps especially, where the inner struc- 
ture is highlighted. Run A is probably more reminiscent of the 
Crab Nebula: the external torus and the two polar jets are well 
visible, as well as the knot (first discovered by Hester et al. 
I1995l in optical Hubble Space Telescope images) and a system 
of wisps, with a brighter arc in what is usually called the in- 
ner ring (Weisskopf et al. 2000). Notice that for this feature we 
find a radius of 0.18 pc, which is reasonably close to the ob- 
served value of 0.14 pc. The torus is associated with the higher 
magnetization equatorial region approximately between 1 and 
3 ly from the pulsar (see again Fig. |3, in particular the emis- 
sion is higher at ~ 2 ly, where the flow converges advecting 
emitting particles coming from higher latitudes (see the corre- 
sponding zone with higher values of 6oo). A feature common to 
both structures is the fading of the emission toward the borders, 
far from the polar axis, where \B' x n'\ becomes small in the co- 



efficient of Eq. fT8t . due to the hypothesis of a purely toroidal 
field (see the discussion below). 

The jets are hollow, since the toroidal field vanishes on the 
axis, but in run A their emission is significant because emitting 
particles are efficiently transported there by the high speed flow. 
This last effect is reduced in run B, where the X-ray jets basi- 
cally disappear, while the diffuse X-ray emission is enhanced 
in the nebula. The explanation is again found by looking at 
Fig. |2| in run B, the hoop stresses are less efficient, due to the 
wider low magnetization region, and particles lose more energy 
through synchrotron radiation before flowing toward the axis. 
It is important to notice that, in both cases, jets are visible in the 
optical band, and actually, at lower frequencies, the difference 
between their appearance in the two runs is reduced. However 
they are not so bright now with respect to their surroundings 
and the diffuse emission from the outer nebula hides them, pre- 
cisely as it happens in the observations. In the central region, 
other than the (overly) bright knot, two boosted arcs are visible 
in run B, reminiscent of the situation encountered in the Vela 
PWN (Helfand et al. l200Tl Pavlov et al. l200Tt . 

Summarizing, it appears that run A, which has a more mag- 
netized wind near the equator, leads to a brighter X-ray torus 
and to a single inner Doppler boosted arc, whereas run B, where 
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run A: surface brightness l„ (optical) run A: surface brightness l„ (X-ray) 
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Fig. 3. Maps of surface brightness 7 y , in logarithmic scale and normalized to the maximum value, for run A (upper row) and 
run B (lower row). Optical images (5364 A) are displayed in left panels, X-ray images (1 keV) on the right. The PWN axis of 
symmetry has been inclined by 30° with respect with the plane of the sky and of 48° with respect to the North, to compare with 
images of the Crab Nebula. For a distance d ^ 2 kpc, we have 1 ly^ 32". 



the low magnetization region around the equator is more ex- 
tended, implies a more diffuse nebular emission and a more 
complex system of rings and arcs in the central region. 

Before concluding this sub-section, let us comment on the 
consequences of the assumption of a purely toroidal magnetic 
field. Although evaluation of the dynamical effects of a poloidal 
field component would require a full 3-D treatment, that is be- 
yond our present goal, the effects of an isotropic field on the 
emission properties are straightforward to check. If we thus 
substitute the term \B' x n'\ in the emissivity with its average 
over the solid angle (compare with Eq. (l23l): 



B' ± = \B'xn'\ = 




we obtain the results shown in Fig. 01 where only case A has 
been considered (similar differences arise for case B). The main 
difference, at a given frequency, between the corresponding im- 
ages in Fig. [21 and |4| is the relative brightness of some of the 



small scale features mentioned above. First of all the knot ap- 
pears to be less dominant than before, as well as the main arc in 
the inner ring. Moreover, as one would expect, the outer parts 
of all ring-like structures are enhanced. In particular, looking 
at the X-ray map, we notice that what we here have referred to 
as the inner ring is now more uniform, as in the observations 
rWeisskoof etal. BDDa . 

4.3. Polarization maps 

Let us now concentrate on the (linear) polarization properties 
of the nebular emission. A preliminary study was reported in 
an earlier paper (Bucciantini et al. 2005b), where maps of the 
polarization fraction (II y ) and direction based on a similar sim- 
ulation were presented. In Fig.|5lwe show the same maps for 
runs A and B, where the central region is zoomed and not 
rotated in the plane of the sky for ease of interpretation (the 
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Fig. 4. Maps of surface brightness for run A corresponding to those in Fig. [3] obtained with an isotropic magnetic field. 



tilt is retained). The two cases are similar: the polarized frac- 
tion is high along the polar axis, where the projected magnetic 
field is always orthogonal to the line of sight, while depolar- 
ization occurs in the outer regions, where contributions from 
projected fields with opposite signs sum up along the line of 
sight. On the right hand side we show the polarization direc- 
tion, with ticks proportional to the fraction displayed on the 
left, together with the corresponding map of surface brightness, 
with the aim of making clear the association between the po- 
larization behaviour and the emission main features. As we can 
see, polarization ticks are basically always orthogonal to the 
toroidal field, displaying the behaviour expected given the in- 
clination of the symmetry axis with respect to the plane of the 
sky. However, in the rings where the velocity is relativistic, de- 
viations of the vector direction due to polarization angle swing 
are also visible (e.g. in the arc of the torus, that is the second 
brightest arc above the central knot). 

Concerning differences between the two runs, we notice 
that the structure of polarization fraction and direction becomes 
more complex for run B, for which also the magnetization map 
in Fig. |2| was more complicated, leading to the presence of mul- 
tiple rings. Notice also that in these zoomed maps it is possi- 
ble to see a small displacement (to the South) of the knot with 
respect to the central position, as in the figures by Hester et 



al. ([1995). Unfortunately, high resolution optical polarization 
maps of the inner region of the Crab Nebula are not available 
yet: all papers mentioned in the Introduction refer to the whole 
nebula, where contributions from the confining medium may 
be significant. Finally, for X-ray polarimetry, which could re- 
ally provide crucial clues to the magnetic structure in the inner 
region, the necessary technology is still to come. 



4.4. Spectral index maps and integrated spectra 

Consider now the spectral properties, starting with the maps of 
spectral index a v in Fig.0 both in the optical and X-ray bands. 
Optical images are obtained by comparing specific intensities 



integrated spectrum 



0.1000 r 




10 19 



Fig. 7. The integrated spectrum as a function of frequency for 
the two runs. Spectra are normalized to X - 5364 A. 

at X - 5364 A and X - 9241 A, wavelengths that have been 
chosen so as to allow direct comparison with the observational 
data presented in the paper by Veron-Cetty & Woltjer ( 1993). 
The X-ray spectral index maps, instead, are obtained by using 
energies of 0.5 and 8 keV, as in Mori et al. ( 120041) . 

The optical images do not show significant variations from 
the injection value of a = 0.6 in the inner region. Once again 
run A seems to be more closely representative of the situation 
in the Crab Nebula, with spectral softening occurring, starting 
from a distance of about 0.5 pc from the central pulsar position, 
in the vicinity of the polar axis, and more gradually at larger 
cylindrical radii. This is consistent with observations, at least 
qualitatively, since softening, in the optical band, is present 
only beyond the torus. A general consideration is that in our 
runs, and especially in run B, the nebular magnetic field in the 
central region seems to be weaker than the value usually as- 
sumed for the average field in the Crab Nebula ~ 3 x 10" 4 G. In 
the torus, our simulations give a field strength lower by a factor 
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run A: polarization fraction n„ (optical) run A: polarization direction + surface brightness (optical) 
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run B: polarization fraction n„ (optical) run B: polarization direction + surface brightness (optical) 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 

Y Y 




0.0 0.2 0.4 0.6 0.8 1.0 -2.0 -1.5 -1.0 -0.5 0.0 



Fig. 5. Maps of (optical) polarization fraction n v (left) and direction (right, superimposed to the surface brightness map), for 
run A (upper row) and run B (lower row). Polarization fraction is normalized against -^j^ 70%, and the ticks length used for 
polarization direction is proportional to the normalized fraction. 



of 2, for run A, and 3, for run B. In addition to that, when com- 
paring to real images, one should also consider that dimensions 
must be rescaled somehow. While in run A the positions of the 
main features are roughly in agreement with the observations 
(see the discussion in Sect. 14.2b . in run B the radii of the in- 
ner ring and the torus are approximately two times larger than 
observed. This naturally concurs to the spectral flatness of the 
corresponding map. 

In the X-rays synchrotron burn-off occurs on much shorter 
distances and the maps are much more shaped: spectral soft- 
ening is much stronger and, as a consequence, basically only 
the inner structure is visible in the adopted range. Notice how 
the values are similar to those measured by Mori et al. ( 120041) 
(where photon index is used for the scale, that is a v + 1), and 
in run A we even observe the brighter (harder) features in the 
polar jets, especially the southern one. 

Finally, let us calculate the integrated spectra (net fluxes) 
as in Eq. (l33l . for the two runs A and B, and compare them 



with spectro-photometric data of the Crab Nebula, at optical 
and X-ray frequencies. In Fig. the net flux F v (in Jansky) is 
plotted against frequency v (in Hertz). Solid lines are the syn- 
thetic integrated spectra obtained from simulations. The four 
crosses in the optical range correspond to data at X - 9241, 
6450, 5364, and 3808 A, corrected for interstellar absorption 
and thermal contribution from the filaments (Veron-Cetty & 
Woltier [T993L whereas the X-ray data point at 1 keV is taken 
from Willingale et al. (2001) and refers to XMM-Newton ob- 
servations. We can notice the different behaviour of the two 
runs. For run A the spectrum steepens because of synchrotron 
burn-off already in the optical band, where the slope is 0.75 
rather than the injection value of 0.6. A spectral break is seen 
at Vb\ « 2 - 3 x 10 15 Hz, as expected (e.g. Veron-Cetty & 
Woltier [T993l> . with the slope changing to 1.15, yielding an in- 
crease of 0.55, comparable to the value 0.5 usually adopted. On 
the other hand, in run B, the lower magnetic field causes losses 
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Fig. 6. Maps of spectral index a v , as defined in Eq. d32l> . for run A (upper row) and run B (lower row). Optical images are 
obtained by comparing specific intensities at X - 5364 Aand X - 9241 Aand are displayed in left panels. X-ray images are 
obtained by using energies of 0.5 and 8 keV, and are displayed on the right. 



to be negligible in the optical and the break to move toward 
X-ray photon energies. 

The most puzzling feature is the presence in both cases of 
further, unexpected, changes of slope, with the appearance of 
what we term as inverse breaks, i.e. flattenings of the spectrum. 
Focusing on the results relative to run A, we see that the spec- 
tral index has a second change from 1.15 to 0.8 at around v^i ~ 
3 x 10 16 Hz; then it increases again to 1.1 at ~ 10 18 Hz; and 
finally decreases to 0.85 at ~ 3 x 10 18 Hz. The behavior of 
run B is qualitatively analogous, although there are differences 
in the position and entity of the changes of slope. 

The reason for this behavior is better understood by look- 
ing at Fig. [8| where we show the cut-off frequency for run A as 
a function of space. In particular, if one looks at the contours, 
it is clear that while the initial break at v&i is due to loss of 
emission from high latitudes, that particles with energies corre- 
sponding to emission at v > vm are not able to reach (consider 
the contours referring to 10 15 , solid line, and 10 16 Hz, dotted 



line, in the figure), the inverse break at v^i appears because 
there is a range of frequencies for which the emitting region 
keeps constant in size. This is strictly related to the particu- 
lar configuration of the magnetic field (see the magnetization 
map in Fig.|3 that is mainly confined to the equatorial region 
due to the presence of large scale vortices. This is confirmed 
by the spatial ditribution of the maximum particle energy 6oo, 
which for high frequencies (thus in X-rays) appears to be con- 
centrated just around the termination shock, so that a further 
increase in frequency does not lead to any further reduction of 
the emitting region. 

The result of the deviations from the expected monotonic 
spectral steepening described above is that X-ray emission is 
too high, by a factor of 2 in run A and 3 in run B. This dif- 
ference of values between the two runs mostly depends on the 
different strength of the magnetic field which causes stronger 
synchrotron losses for run A than for run B, although still not 
enough to explain the observations quantitatively. 
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run A: maximum particle frequency 
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Fig. 8. The distribution of the maximum particle frequency (in 
Hz units), that is the critical frequency of Eq. (TT9t calculated 
ignoring line of sight effects. Contours refer to the values of 
10 15 Hz (solid line), 10 16 Hz (dotted line), 10 17 Hz (dashed 
line), 10 18 Hz (dot-dashed line). 



5. Conclusions 

In this paper a complete suite of diagnostic tools for deriv- 
ing synchrotron emission properties from numerical relativis- 
tic MHD simulations has been presented. The method is quite 
simple and general, so that it can be applied to any numeri- 
cal scheme and to any astrophysical source of synchrotron ra- 
diation. The application presented here concerns the emission 
from PWNe, based on the axisymmetric model by Del Zanna 
et al. d2004l» . In this case, synchrotron radiation is the primary 
diagnostic for the physics in the source and may give impor- 
tant clues about the structure of the pulsar wind itself, since 
its properties are found to be reflected in the nebular emis- 
sion. In particular we have studied the differences arising from 
two different magnetic structures of the wind: one in which the 
striped (low magnetization) wind region is very narrow, and 
the toroidal field has a maximum near the equator (run A), and 
another where the maximum occurs at a much higher latitude 
above the equator (run B). In both runs the average magnetiza- 
tion was the same, <x e ff ~ 0.02, and high enough to induce the 
presence of supersonic polar jets. 

The width of the striped wind region of low magnetization 
is usually considered as directly linked to the inclination an- 
gle between the pulsar's rotation and magneto spheric axes, thus 
such a study could yield crucial information otherwise difficult 
to obtain. By looking at surface brightness maps, especially in 
X-ray, it appears that run A produces images which are re- 
markably similar to that of the Crab Nebula observed by the 
ROSAT (Hester et al. 113951 . Chand ra (We isskopf et al.[2000), 
and XMM-Newton (Willingale et al. 2001 ) satellites. An equa- 



torial torus and two polar jets are clearly visible, together with 
an inner ring, with Doppler boosted emission concentrated in 
a main arc, and a bright knot just SE from the pulsar position. 
Particularly noticeable is the presence of the jets, necessarily 
hollow because of our assumption of axisymmetry and of a 
purely toroidal field (though the case of an isotropic field has 
been treated here too): these did not appear as emitting features 
in previous simulations. Run B produces instead maps which 
resemble vaguely the Vela PWN (Helfand et al. 200 lj Pavlov et 
al. l2Q0Hl2Q03l> . with a couple of bright arcs NW from the cen- 
tral knot in the inner region. This is entirely due to the different 
flow structure determined by the lower magnetization around 
the equator (hoop stresses and small scale vortices develop at 
higher latitudes). 

In the present paper we have used observations mainly of 
the Crab Nebula as a benchmark for our diagnostic techniques, 
for the simple reason that the Crab Nebula is by far the best 
studied PWN. We have briefly mentioned the Vela case too, 
but one may wonder whether by tuning our model parameters 
the resulting surface brighness maps could reproduce also dif- 
ferent features encountered there or in some other PWN. One 
of the most puzzling problems has always been the jet in the 
Vela PWN, which is much brigher on the wrong side of the 
pulsar position (that for which the jet is receding from us). The 
best explanation is probably to invoke kink-type instabilities 
that manage to bend sections of the counter-jet toward us, thus 
enhancing their emission, but in order to test this idea full 3- 
D simulations would be required. Regarding other objects, we 
may consider the case of the PWN around pulsar B 1509-58. 
In the central region multiple knots are observed (Gaensler et 
al. 120021) and the brightest ones are located between the pulsar 
position and the boosted arc on the torus, opposite to the direc- 
tion of the (main) jet. Again, these features cannot be explained 
within our 2-D model, possibly they could arise due to time 
dependent small-scale Kelvin-Helmholtz instabilities, as it has 
been suggested for the variable optical wisps and sprites in the 
Crab Nebula (Begelman [T999l Bucciantini & Del Zanna 2006 ), 
but in any case the approximation of axisymmetry must be re- 
laxed. In the same source, the ring-like region named RCW 89 
could be produced by the counter-jet (very faint in the Chandra 
images) hitting the contact discontinuity between the PWN and 
the expanding ejecta. This scenario is well within our assump- 
tions and its accurate modeling is left as future work. 

We have also computed detailed polarization maps. This 
is expected to be a very powerful diagnostic technique for the 
magnetic structure of the inner regions of PWNe, since syn- 
thetic surface brightness maps are likely to be too largely dom- 
inated by Doppler boosted features, in the vicinities of the cen- 
tral source, to be sensitive to small changes in the magnetic 
field strength and direction (see Bucciantini et al. 2005b] for a 
wider discussion). The results show a high degree of linear po- 
larization close to the axis of symmetry (not far from the theo- 
retical maximum of 70%, given the spectral index a = 0.6), 
and a strong depolarization in the outer equatorial regions. 
Relativistic effects related to Doppler boosting (polarization 
angle swing) are also visible in the brighter arcs. High resolu- 
tion optical, not to mention X-rays, polarization observations, 
needed for a comparison of the numerical results with the real 
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data, aimed at probing the inner PWN magnetic structure, are 
lacking, at present, but are highly desirable given their diagnos- 
tic potential. 

Synthetic spectral index maps have been computed here for 
the first time on the basis of MHD numerical simulations. At 
a qualitative level, the results of our run A are in good agree- 
ment with both optical (Veron-Cetty & Woltier [T993l) and X- 
ray (Mori et al. 2004 ) data: spectral softening is weak in the 
central region and stronger at the borders of the torus, espe- 
cially in the X-rays, where even the harder features inside the 
jets highlighted in Chandra data are reproduced. However, the 
nebular field appears to be too low to explain the observations 
quantitatively, giving too weak a softening, especially for run 
B. This is further confirmed when looking at the integrated 
spectra, which show that the spectral break is too weak and 
occurs at a too high frequency. The result is that, if we nor- 
malize the spectra to a flux in the optical, X-ray data are over- 
estimated by a factor 2 for run A and 3 for run B. We have 
attributed this behavior to the distribution of the nebular mag- 
netic field and, consequently, to that of the maximum energy of 
emitting particles which serves as a local synchrotron cut-off. 
Such distributions basically show that the magnetic field is too 
strongly affected by the compression due to large scale flow 
vortices and piles up around the termination shock. We have 
also verified that this situation invariably occurs for different 
values of the simulation parameters. 

In spite of these successes in PWN modeling, in order to 
proceed further 3-D simulations are certainly required. The 
strongest limitation appears to be the assumption of a purely 
toroidal magnetic field, both for the overall dynamics (e.g. al- 
lowing for kink-type instabilities along the axis and possibly 
reconnection at the equator) and for the emission properties, 
where a certain amount of disordered magnetic field on small 
scales could help to reproduce the observations, as we have 
shown in maps for the optical and X-ray bands in the case of an 
isotropic field (see also the discussion in Shibata et al. (2003). 
Another issue not properly addressed here is time variability: 
though the overall evolution of the PWN / SNR system is nearly 
self-similar in the free expansion phase, the small scale vortices 
around the termination shock are variable on a time scale of a 
few years (see also Bogovalov et al. 2005). We leave the dis- 
cussion of these two important issues to future papers. 
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